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The orientational dependence of charge carrier mobihties in organic semiconductor crystals and the 
correlation with the crystal structure are investigated by means of quantum chemical first principles 
calculations combined with a model using hopping rates from Marcus theory. A master equation 
approach is presented which is numerically more efficient than the Monte Carlo method frequently 
applied in this context. Furthermore, it is shown that the widely used approach to calculate the 
mobility via the diffusion constant along with rate equations is not appropriate in many important 
cases. The calculations are compared with experimental data, showing good qualitative agreement 
for pentacene and rubrene. In addition, charge transport properties of core-fluorinated perylene 
bisimides are investigated. 



I. INTRODUCTION 

Due to their low production costs and easy processabil- 
ity organic semiconductor devices are promising materi- 
als for organic light emitting diodes (OLEDs)^' — organic 
field effect transistors (OFETs),''"^ radio frequency iden- 
tification tags (RFIDs)^'^ and solar cells, ^"^"^ to mention 
just a few. The performance of these devices depends cru- 
cially on the charge transport. Therefore, it is important 
to understand the basic principles of charge transport in 
these materials. 

Various models have been proposed which are often 
contradictory. The band theory, which is well estab- 
lished for inorganic covalently bonded materials, is not 
particularly appropriate for organic conductors, because 
organic molecular crystals are only weakly bound by van 
der Waals interactions causing the molecules to be much 
more flexible. Due to the complex nodal structure of 
the molecular orbitals the transfer integrals between the 
monomers are very sensitive to even small nuclear dis- 
placements. That is why lattice vibrations play a more 
important role in organic than in inorganic materials, as 
they destroy the long range order and lead to a charge 
carrier localization.^"* To account for these vibrations, 
a variety of models have been proposed which incorpo- 
rate the local (Holstein)^^ and the nonlocal (Peierls)^^ 
coupling. The latter leads to a polaron model where 
the charge carrier is partially localized and dressed by 
phonons.*^"^- The fluctuations of the coupling between 
the molecules are of the same order of magnitude as 
the average coupling, leading to a rather strong local- 
ization. Other models have been suggested, where the 
charges are assumed to be localized and the inter- and 
intramolecular vibrations are treated classically<^"— 

At higher temperatures, it is often appropriate to as- 
sume that the charge is localized due to the thermal disor- 
der of the molecules and that charge transport occurs via 
thermally activated hopping. In some cases room tem- 
perature should be sufficient for this assumption to be 
justified. We apply this hopping model to study the de- 



pendence of the charge carrier mobility on the molecular 
structure and morphology as well as its angular depen- 
dency. The latter point is important since most organic 
crystals show a pronounced anisotropy for the transport 
parameters which has to be taken into account for device 
design. Furthermore, it is known that the mobility is 
very sensitive to the arrangement of the monomers and 
that already small changes in their alignment can alter 
the transport parameters dramatically/^ 



A promising class of materials for organic electronics 
are perylene bisimides. Due to their light resistance^! 
and intense photoluminescence^^ they are widely used 
as robust organic dyes in the automobile industryi^l Fur- 
thermore, they show a considerable electron mobility22ri^ 
and a high electron affinity, '^^-'^ That is why they serve 
as n-type semiconductors for organic field effect transis- 
tors^"— and as electron acceptor material in organic so- 
lar cells,^-ia 



Section |TT] describes the theoretical background of the 
applied model as well as details of the numerical calcula- 
tions and computational approaches. It is shown that the 
master equation approach is particularly faster than the 
well-known Monte Carlo method. Furthermore we eluci- 
date why the commonly applied approach to calculate the 
mobility via the diffusion constant along with rate equa- 
tions^ is not appropriate in many important cases. In 
Sec. nil Al we consider the frequently disputed question 
if the Einstein relation holds even for more disordered 
(amorphous) materials.—"— In Sec. IIII Bl we show re- 
sults for the orientational and morphological dependency 
of the mobility for pentacene, rubrene and two fluori- 
nated perylene bisimides. The first two materials are ex- 
perimentally and theoretically well investigated^'^'^"— 
which allows for the comparison with experimental data. 
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II. THEORY AND MODELING 

A. The Marcus hopping model 

In this work, a hopping mechanism is assumed for the 
motion of the charge carriers. The hopping rate from a 
site i to j is given by the Marcus equation^i^ 



XksT 



exp 



4AfcBT 



(1) 



where Vji is the electronic couphng parameter, A is the 
reorganization energy, T is the temperature, is the 
Boltzmann constant and h = /i/(27r) where h is the 
Planck constant. The energy difference AEji between 
the two hopping sites is caused by an external electric 
field F. If the material is less ordered or even amorphous, 
each molecule experiences slightly different surrounding 
effects (such as polarization) that lead to different site 
energies i??. These energy differences furthermore con- 
tribute to AE 



AE,, = {E^ E^) - qF, 



(2) 



where q is the charge which equals the positive or neg- 
ative unit charge and r*,; is the distance vector between 
sites i and j. Marcus rates have been used before for cal- 
culating the anisotropy of the charge carrier mobility?^ 
but with AEji = 0. 

The interaction of the charge carriers with the phonons 
is partially considered by the reorganization energy. Due 
to the weak van der Waals interactions between or- 
ganic molecules, it can be divided into an internal (in- 
tramolecular) and an external (intermolecular) part, i.e. 
A = Aint + Aext- The intramolecular reorganization en- 
ergy Aint is due to the geometry changes of the donor and 
the acceptor monomer upon the charge transfer process. 
The external reorganization energy Aext covers the ener- 
getic changes concerning the surrounding, caused by lat- 
tice distortion and polarization. For oligoacenes Aext was 
shown to be about one order of magnitude smaller than 
Aint-^^'^ Furthermore, is was demonstrated that Aint of a 
molecule is lower in a cluster than in gas phase and that 
the total reorganization energy of naphthalene is closer 
to Aint in the gas phase than to Aint in the cluster That 
is why the external reorganization energy is neglected in 
this paper and the internal reorganization energy of the 
monomer in vacuum is used for A. 

The Marcus theory was originally derived for outer 
sphere electron transfer in solvents. It stems from 
time dependent perturbation theory (Fermi's Golden 
rule) and describes a non-adiabatic charge transfer where 
the charge carrier is localized at the donor or acceptor 
molecule respectively. Treating the coupling as a pertur- 
bation requires that Vji is small compared to A/4, which 
corresponds to the activation energy for the charge car- 
rier to change place (for AEji = 0). Furthermore, the 
thermal relaxation (the geometric reorganization) has to 



be fast in comparison with the transfer so that the system 
can be assumed to be in thermal equilibrium during the 
transfer. In addition, the theory is restricted to the high 
temperature case since tunneling is neglected completely 
and the molecular vibrations are treated classically, what 
requires k^T ^ fko. These restrictions of the Marcus 
theory in the context of charge transport are discussed 
elsewhere i^i^ Despite all imperfections it is widely used 
for charge transfer in organic crystals^"— i^"— and one 
can certainly assume that this theory is suitable for the 
purpose of a qualitative charge transport analysis. 



B. The master equation approach 

The master equation approach was used to describe 
the transport process. In the case of low charge carrier 
densities, the master equation, which describes the hop- 
ping of the charge carriers in the organic semiconductor, 
has the simple linear formal 



dpi 
dt 



(3) 



where pi denotes the probability that the lattice site i 
is occupied by a charge carrier. The index j sums over 
all other sites. In principle, it is also possible to include 
repulsive forces between the charge carriers in the master 
equation in order to account for higher charge carrier 
densities. However, in the case of low densities, even the 
quite simple Eq. ^ leads to good results. 

In the steady state, a dynamic balance is reached where 
the occupation probabilities for the sites do not change 
anymore and dpi/dt in Eq. ([3]) equals zero. Since this 
equation holds for all sites in the crystal, this results in 
a linear system of equations. 



N -p = 0. 



(4) 



p contains the unknown pi and N is a negative semidef- 
inite sparse matrix that contains all hopping rates Vji. 
For one dimension N is 
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(5) 

The columns correspond to the initial sites i of the charge 
carrier and the lines correspond to the final sites j, i.e., 
the jump rate Vji from i to j appears in the iih. column 
and the jth line. The diagonal elements contain the neg- 
ative sum of all hopping rates away from the respective 
site. 

The infinite matrix N is approximated by a finite ma- 
trix with cyclic boundary conditions, i.e., a charge carrier 
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that leaves the crystal at one side reenters at the oppo- 
site side. This means for the example matrix depicted in 
Eq. ([U that the charge which jumps from site 4 in posi- 
tive direction ends at site 1. For this boundary condition 
to be applicable it has to be assured that the hopping 
rate from site 4 to site 1 in negative direction is negligi- 
ble. This results in a constraint for the minimum size of 
the matrix. 

The matrix in Eq. ([5]) was extended to three dimen- 
sions resulting in a (Sn^rim) x (Srtrfrtm) matrix where rt^ 
is the number of unit cells in each direction and rim is 
the number of monomers per unit cell. In this work all 
monomers within a cube of three unit cells length in each 
dimension of the crystal are taken into account. It was 
verified that a bigger matrix with more than 3x3x3 
unit cells does not change the result. The hopping rates 
were calculated from one monomer to all other monomers 
in the same and in the adjacent cells. Since the jump 
rate, Eq. ([1]), implicitly depends on the distance via the 
electronic coupling Vji, larger jump distances can be ne- 
glected. 

Solving Eq. ^ and taking into account the normaliza- 
tion condition J^i Pi = 1 provides the occupation proba- 
bilities for all sites. (For AEji = 0, it is the same for all 
sites.) These probabilities can then be used to calculate 
the mobility of the charge carriers in field direction from 



F 



with the average velocity 



(6) 



(7) 



where Vi is the resulting velocity at site i, 
is the average displacement at site i in field direction and 



(8) 



(9) 



is the dwell time of the charge carrier at site i. Equa- 
tions dH) to dll) result idJ^ 




12 j ^ji {^ji 
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F 



(10) 



In order to simplify the calculation of the mobility 
within such a jump rate approach, the mobility is often 
calculated without external field because the occupation 



probabilities of the sites do not differ in this case and one 
does not have to solve the master equation (U). Since 
Eq. (fTO|) is not applicable in that case (because F = 0), 
the mobility is calculated via the diffusion coefficient D 
and the Einstein relation^'^ 



-D. 



(11) 



Different equations are found in the literature^^— 
to evaluate D. Considerations similar to those above for 
the mobility seem to provide 



2ndV ' 



2n E^ 



(12) 



where n is the spatial dimensionality. Since the diffusion 
is regarded in one dimension here, n equals 1 and 



1 



D = - 



2^' 



where 



^3 



(13) 



(14) 



is the variance of the charge carrier position at site i in 
the direction of the unit vector e. Equations (|9]) and (fT2)) 
to (11411 finally result in 



(15) 



for the diffusion coefficient in the direction of e. It is 
worth mentioning that Eq. ([T5|) holds even in the pres- 
ence of an external field (see Appendix |X| . 

Without external field and assuming that all lattice 
sites are equal (i.e. AEji — 0), the last equation simplifies 



to 



44.75.76 



(16) 



It is important to note, that the diffusion constants in 
Eqs. ([T5|) and are not strictly correct. Just if the unit 
cell of the crystal contains only a single molecule and if 
the crystal structure is perfectly translation-symmetric, 
i.e. E'^ = Ej for all monomer pairs, cf. Eq. ([2]), these 
equations become correct. 

However, in less ordered or even amorphous materi- 
als the site energies E^ and Ej are different because of 
the differing surroundings for each lattice site. In that 
case, the occupation probabilities pi differ and the mas- 
ter equation has to be applied. In the case of strongly 
different Ef, even Eq. (|15p becomes incorrect since the 
charge carrier can be "trapped" between two lattice sites 
with similar energyj^^ see Fig.[T^: Because of the energet- 
ically unfavorable surrounding, the charge carrier jumps 
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FIG. 1. The charge carrier is "trapped" between two lattice 
sites, a) The surrounding of the monomer causes an energetic 
"pit", b) Strongly differing jump rates lead to a capturing. 



back and forth between the same sites all the time. These 
moves do not contribute to a macroscopic spreading of 
the occupation probability of the charge carrier with the 
time. That is why the averaging in Eq. (jlSp overesti- 
mates the true macroscopic diffusion coefficient. This 
problem does not appear in Eq. (|10p since the fji is not 
squared as in Eq. (fT5|) . For that reason the contribution 
of the trapped charge cancels when summing over all lat- 
tice sites. And even in perfectly ordered crystals where all 



jump rates are symmetric, i.e. Vji 



(without external 



field) , such a trapping can occur if different sites exist in 
the elementary cell of the crystal and if the hopping rates 
within the cells differ from those to neighbored unit cells, 
see Fig.[T|D: Here, the charge carrier jumps back and forth 
between two monomers with a high coupling because the 
coupling to the other neighbors is lower. In such cases 
Eq. (fTO]) in conjunction with Eq. ([TT|) provides correct dif- 
fusion coefficients while Eqs. ([TSj) and ([T6|) overestimate 
the values for D. 



C. The Monte Carlo approach 

The master equation results were verified with Monte 
Carlo simulations applying the algorithm of Houili et 
al.j^ but without any interaction between the charge 
carriers. The mobility and the diffusion coefficient were 
calculated via 




and 



D = ((r,je 



(17) 



(18) 



respectively. The time dependent average position {rjij^) 
and the variance {{rjie— {r.jie))^) have been averaged 
over a sufficient number of simulation runs to obtain 
smooth lines. It was checked that both average and vari- 
ance show a linear time dependence in order to secure 
the stationary state. 



The Monte Carlo approach is just an alternative way 
to solve the master equation ([3]). It is a feasible way 
to log the atomic scale motions underlying the trans- 
port properties as a function of time. However, as this 
is a stochastic method, many simulation runs are needed 
in order to achieve an acceptably low statistical error 
such that sufficiently significant values are obtained for 
the mobility and the diffusion coefficient. Furthermore, 
one has to take care that the stationary state is reached 
within the simulation time. This is a serious problem in 
the case of strongly disordered materials. In contrast to 
that, the approach used here by solving the matrix equa- 
tion ([4]) which provides the stationary state by means of 
analytic numerical methods guarantees the stationary so- 
lution and is furthermore numerically more efficient than 
Monte Carlo simulationsi^ 



D. Quantum chemical methods 

The electronic coupling Vji and the reorganization en- 
ergy A needed for the hopping rate, Eq. ([1]), are deter- 
mined by quantum chemical first principles calculations. 
In order to calculate A, the geometry of the isolated 
monomer was optimized for the charged and the neu- 
tral state. The energies Eq and Ec of the neutral and 
the charged monomers in their lowest energy geometries 
and the energies E^ and E* of the neutral monomer with 
the ion geometry and the charged monomer with the ge- 
ometry of the neutral state are calculated to get the in- 
tramolecular reorganization energy^ 



A = A, + Ao = [El - E,) + {E* - Eo) 



(19) 



cf. Fig. [2] For all quantum chemical calculations the 
TURBOMOLE program package^^ was used. The cal- 
culations were conducted via density functional the- 
ory using the hybrid generalized gradient functional 
B3-LYP— with the correlation consistent polarized 
valence double zeta basis set (cc-pVDZ)2^ for all 
atoms. This functional was chosen because it has been 
shown that it leads to quite good results for describ- 
ing the ionization-induced geometry modifications of 
oligoacenes 

The electronic couplings were calculated as described 
by Li et al^ resulting in 



^ji ~ ^{Hji + Hjj)Sj. 



(20) 



with 



{ipj\HKs\Vi) 



For hole (electron) transport ipi and ipj are the HOMO 
(LUMO) orbitals of the respective isolated monomers and 
Hxs is the Kohn-Sham operator of the neutral dimer 
system. Ha and Hjj are the site energies of the two 
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FIG. 2. The potential energy surfaces of the neutral and the 
charged monomer. The dashed arrows indicate the vertical 
transitions from one state to the other. Aq and Ac are the two 
contributions to the reorganization energy, see Eq. p9|) . 



The first exponential function describes the tunneling of 
the charge and the Boltzmann-type exponential function 
accounts for thermally activated jumps upwards in en- 
ergy. Hops to lower energies are not thermally activated. 

A simple cubic lattice of sites with a lattice constant of 
1 nm was used. In order to achieve a sufficient statistics 
for the site energies the lattice consisted of 80 x 40 x 40 
sites. For a given site only the hops from and to the 
26 adjacent sites were considered. Calculations with a 
bigger lattice and also further jump targets taken into 
account did not affect the result. 



III. RESULTS AND DISCUSSION 
A. Validity of the Einstein relation 



monomers, Sji is the spatial overlap and Hji is the charge 
transfer integral in the non-orthogonalized basis. 

The arrangement of the monomers in the crystal was 
extracted from X-ray crystal structure data which was 
retrieved from the Cambridge Structural Database. 

E. The Gaussian disorder model 

It has been argued that the Einstein relation, Eq. (ITT]) , 
does not hold in disordered organic materials in gen- 
eral^"— or at least if additionally an external field is 
applied In fact it turned out that this is only true 
for rather high charge carrier densities;^ low tempera- 
tures and high electric fields which are out of the scope 
of the present work. At extremely low temperatures, the 
thermal energy of the charge carriers is not sufficient to 
reach sites which are higher in energy and only energy- 
loss jumps occur. In that case, neither nor D depends 
on the temperature.!^ For low fields, the transport coef- 
ficients are independent of the field, ^^'^'^ but for higher 
fields nonlinear effects become important and D/ix in- 
creases with increasing fieldf^i 

A strongly disordered organic semiconductor was sim- 
ulated by means of the Gaussian disorder model^^ with 
a Gaussian shaped density of states, 

where the standard deviation a is called the energetic 
disorder of the simulated material, in conjunction with 
the Miller- Abrahams jump rate^ 

Vji = iyoexp(-27rji) x (22) 

/exp(-^), A£;,,>o 
1 1, A^j, < 

where = 10^"^ s~^ is the attempt-to-jump frequency 
and 7 = 5- 10^ m^^ is the inverse localization radius. 



The mobility and the diffusion coefficient were calcu- 
lated by the master equation approach in conjunction 
with the Eqs. PH]) and and by the Monte Carlo ap- 
proach using Eqs. (flTl) and (flSl) respectively. The Gaus- 
sian disorder model described in Sec. Ill El was used. In 
the Monte Carlo simulation, the average and the vari- 
ance of the charge carrier position has been averaged over 
50.000 trajectories and the simulation time has been up 
to 1 s. 

Figure|3]shows the results as a function of the energetic 
disorder cr, cf. Eq. (PTT) . The mobility varies over several 
orders of magnitude and the results of Eqs. ([TU]) and p7)) 
match exactly. This is not the case for the diffusion coef- 
ficient calculated with Eq. ([T5|) and (ITSl) . With increasing 
energetic disorder, the deviations between these two ap- 
proaches to calculate D increase. These deviations are 
not caused by the field because for cr = the results 
match. In order to decide which one is the right ap- 
proach, the ratio D/fi is plotted as well. One clearly 
sees that in the case of Monte Carlo the Einstein re- 
lation, Eq. (fTTj) . is valid, whereas D/fi calculated with 
Eqs. (|15p and (llOp deviates from the Einstein relation. 
The two mobility equations lead to the same results. 
Thus, Eq. (ITSj) and also the frequently used Eq. (ITB|) pro- 
vide incorrect diffusion constants for energetically inho- 
mogeneous materials. In any case it is advantageous to 
employ the master equation in conjunction with Eq. (|10p 
to calculate the mobility as this provides correct diffusion 
constants without numerical noise and with low compu- 
tational demands. 



B. Angular dependence of the mobility in crystals 

If not otherwise stated, the calculations have been con- 
ducted with an electric field of 10'' V/m and a tempera- 
ture of 300 K. The molecules under investigation are de- 
picted in Fig. |T| and the crystallographic parameters of 
the corresponding crystals are listed in Tab. HI 
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FIG. 3. (Color online) Mobility, /i, (left) diffusion coefficient, D, (middle) and the ratio D/fj, (right) as a function of the energetic 
disorder, calculated with the rate equations (|IO|l and (|f 5|) respectively, and via Monte Carlo simulation. The calculations were 
conducted at T = 300 K and F = 10^ V/m. 




a) b) c) d) 



FIG. 4. The molecules investigated in this work: a) pen- 
tacene, b) rubrene, c) PBI-F2, d) PBI-(C4Fg)2. 



TABLE I. Lattice constants and angles for the unit cells of 
all calculated crystals. 





a [A] 


b[A] 


c[A] 


a n 


/3 n 


7 n 


Ref. 


pentacene 


6.27 


7.78 


14.53 


76.48 


87.68 


84.68 


50 


rubrene 


26.86 


7.19 


14.43 


90.00 


90.00 


90.00 


58 


PBLF2 


17.46 


5.28 


15.28 


90.00 


110.90 


90.00 


32 


PBL(C4F9)2 


10.57 


12.89 


16.68 


66.86 


76.52 


84.62 


96 



1. Pentacene 

Pentacene (see Fig. exists in several morpholo- 
gies. Here the structure described by Mattheus etali^ 
(at 293 K) was investigated. The unit cell contains two 
differently orientated monomers. Pentacene is known to 
be a hole conductor, but for comparison, the electron 
transport is regarded here as well. The reorganization 
energy was calculated to 92meV for holes and 131 meV 
for electrons. This is in good agreement with values re- 
ported before (98 and 95meV for holes^^^ and 132 meV 
for electrons 4^) 

Figure [5] shows the mobilities of holes and electrons in 
the crystal in all three dimensions. For better legibility 
Fig. [6] shows two dimensional cross sections orthogonal 
to the a*, b* and c* direction respectively. The magni- 



TABLE II. The most important electronic couplings and the 
reorganization energy in the pentacene crystal for electrons 
and holes, cf. Fig. [T] 





h+ [meV] 


e [meV] 


Vi 


90.69 


85.18 


V2 


55.05 


89.66 


Vs 


39.68 


50.00 


Vi 


36.62 


47.10 


A 


92 


131 



tudes of the hole and electron mobility are quite similar. 
For both types of charge carriers the transport is almost 
two dimensional since the minimal mobility, that is found 
in the c* direction, is very low (0.2cm^/Vs for holes and 
1.3cm^/Vs for electrons) compared with the other direc- 
tions. This can be explained by the electronic couplings. 
The highest ones are listed in Tab. [Ill The directions of 
the corresponding charge transitions are drawn in Fig. [T] 
All of them are coplanar in the ab plane. For holes, the 
biggest coupling belonging to a transition with a compo- 
nent in c direction is one order of magnitude lower than 
the lowest coupling listed in Tab. |lT] (electrons: about 
factor 5 smaller). The highest couplings for holes belong 
to the transitions in [110] direction, the second highest to 
the [110] direction. The reverse is true for electrons. That 
is why the directionality of the mobilities for holes and 
electrons differ in the ab plane. The maximum mobility 
for holes (18.5cm^/Vs) is found at 132°, the maximum 
for electrons (13.7cm2/Vs) at 37°. 

Figure [S] shows a comparison between the calculation 
and some experimental mobility values for holes. ^'^ Please 
note, that the crystal orientation could not be determined 
in the experiment The measured mobility varies be- 
tween 0.66 and 2.3cm^/Vs. This shows that the calcu- 
lated maximal mobility is almost one order of magnitude 
too big. However, in highly purified single crystals of 
pentacene a mobility of 35 cm^/V s has been measuredi^ 
It was also experimentally confirmed that the mobility in 
the ab plane is much larger than along the c* axis;"^^ This 
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FIG. 5. (Color online) The mobility for holes (top) and elec- 
trons (bottom) in the pentacene crystal for F = 10^ V/m and 
T = 300 K. 

is in agreement with our calculations where the minimal 
mobility of about 0.2 cm^/V s is in c* direction. For room 
temperature and lower, the measurements showed a tem- 
perature dependence of the mobility following cx T^" 
with a positive n indicating band transport.^ - While this 
is not in accordance with the thermally activated hop- 
ping model used here, it was also shown that above 
room temperature a different transport mechanism dom- 
inates the mobility. A further reason for the overestima- 
tion of the mobility is that the nonlocal electron-phonon 
couplingiir— is neglected in our model. While the ab- 
solute values do not match the measured mobilities, the 
qualitative dependency on the crystal direction fits to the 
experimental results. 

2. Rubrene 

Rubrene (see Fig. |3}d) is a hole conductor. It crystal- 
lizes with four differently oriented monomers in the unit 
cell. The calculations were conducted using the morphol- 
ogy described by Jurchescu et al^ at 293 K. Table Hill 
shows the reorganization energies and the values of the 



TABLE III. The most important electronic couplings and the 
reorganization energy in the rubrene crystal for holes and elec- 
trons, cf. Fig. [8] For comparison calculated values for holes 
from Refs. [U andH are shown. 





h+ [meV] 


e [meV] 


h+ [meV}^ h+ 


[meVF 


Vi 


95.73 


49.40 


89 


83 




16.38 


5.55 


19 


15 




1.36 


0.59 






Vi 


0.24 


0.24 






X 


146 


199 


152 


159 



four highest electronic couplings. The couplings next in 
size are two orders of magnitude smaller than the smallest 
coupling listed. This is in agreement with previous cal- 
culations 4^*^ The hopping paths corresponding to these 
couplings are drawn in Fig. |8l The largest coupling {Vi) 
is between equally oriented monomers along the b direc- 
tion, which is the smallest lattice constant. The second 
largest couplings are between monomers which lie in the 
same plane perpendicular to the a axis. V3 is the cou- 
pling between these planes and V4 is the coupling between 
monomers in the same plane perpendicular to the b axis. 

In contrast to pentacene, the electronic coupling for 
holes and electrons in rubrene differs remarkably. That 
is why the calculated mobility for electrons is about one 
order of magnitude smaller than for holes, see Fig. HI But 
unlike pentacene, the angular dependence of the mobility 
is qualitatively the same for both types of charge carri- 
ers. For holes a three dimensional depiction is shown in 
Fig. [ini The maximum mobility (20cm^/Vs for holes 
and 3cm^/Vs for electrons) is in b direction because of 
the short lattice constant in that direction and the re- 
sulting strong electronic coupling. The lowest mobility 
(0.03cm^/Vs for holes and 0.003 cm^/Vs for electrons) 
is in a direction. The main contribution to the mobil- 
ity in that direction are the zig-zag jumps between the 
planes perpendicular to b which are marked with V3 in 
Fig. [5] and the zig-zag jumps between the planes perpen- 
dicular to the c axis marked with V4. The corresponding 
couplings are more than one order of magnitude smaller 
than the next highest coupling V2. The zig-zag jumps 
corresponding to V2 are the main contribution to the mo- 
bility in c direction. 

Figure [S] shows some experimental mobility values for 
holes for the ba planei^'^'^ As for pentacene the calcu- 
lation overestimates the mobility. The calculated max- 
imum mobility is four times larger than the measured 
value. The mobilities for pentacene and rubrene calcu- 
lated in Ref. i41| with a similar approach seem to fit better 
to the experiment. Yet it seems that in their calculation 
a wrong dwell time of the charge carriers was used (cf. 
Sec. Ell. 

The reorganization energy for rubrene is much higher 
than for pentacene. It was shown that the low-frequency 
bending of the phenyl side-groups in rubrene around the 
tetracene backbone contributes strongly to A.-^ How- 
ever, this bending might be impeded in the crystal and 
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FIG. 6. (Color online) The mobility for holes and electrons in the pentacene crystal in the ab plane (left), ac plane (middle) 
and be plane (right). The parameters are the same as in Fig. [S] For comparison some experimental values-- are plotted. Note 
that in the experiment the crystal orientation could not be determined^ and therefore the experimental data is rotated to fit 
best. 




FIG. 7. (Color online) The most important hopping paths in 
the pentacene crystal. Direction of view is parallel to the c* 
axis. 



a smaller reorganization energy would lead to an even 
higher mobility. 

Temperature-dependent measurements in rubrene 
have shown a decrease of the mobility around room tcm- 
perature i^^'^^ This is an indication for band transport. 
However, the qualitative anisotropy of the mobility cal- 
culated with the hopping model fits quite well to the 
measurements. 



3. PBI-F2 

The core-fluorinated perylene bisimide PBI-F2 de- 
scribed by Schmidt et ali^ and depicted in Fig. ^ was 
analyzed. This material is quite interesting for appli- 
cation since it is remarkably air stable because of its 
electron-withdrawing substituents which makes the elec- 
trons less susceptible to trapping with oxygen. The pla- 




FIG. 8. (Color online) The most important hopping paths in 
the rubrene crystal. Direction of view is parallel to the a axis 
(left) and the b axis (right) respectively. The black and the 
grey monomers have a different position in b direction. 



narity of the perylene core is only slightly distorted by 
the core fluorination which leads to a torsion angle of 
3°!^ It was shown that PBI-F2 has a narrower valence 
band and a broader conduction band than the unsub- 
stituted PBI, mainly due to the altered molecular pack- 
ing.— The unit cell contains two differently orientated 
monomers. In contrast to pentacene and rubrene, PBI-F2 
is an electron conductor which is caused by its high elec- 
tron affinity. The electronic couplings for electrons and 
holes differ remarkably. The strongest couplings are col- 
lected in Tab. IIVI The couplings which are not listed are 
at least one order of magnitude smaller than the small- 
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FIG. 9. (Color online) The mobility for holes and electrons in the rubrene crystal in the ba plane (left), ac plane (middle) and be 
plane (right). The parameters are F — 10^ V/m, T = 300 K. For comparison some experimental values for hole mobilitie a^'^'^^'^° 
are plotted for the ba plane. 
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FIG. 10. (Color online) The mobility for holes in the rubrene 
crystal in all three dimensions. The parameters are the same 
as in Fig. [9] 



est coupling mentioned. The strongest coupling for elec- 
tron transport is found between monomers shifted along 
the b direction, see Fig. [TT] Note that this is about 300 
times bigger than the coupling next in size, which is the 
one between two differently orientated monomers within 
the same unit cell. The result is an almost one dimen- 
sional charge transport along the b direction, see Fig. [T^] 
and [131 This might be problematic for application, since 
the charge transport gets very sensitive to lattice distor- 
tions, because the electron cannot easily pass at lattice 
defects which cannot be avoided in real crystals. 

Whereas the coupling between b shifted monomers is 
very strong for electrons, this is surprisingly not the case 
for holes. Their coupling is more than two orders of mag- 
nitude smaller than the electron coupling. This is con- 
firmed by other calculations The reason can be found 
in the differing nodal structure of the HOMO and the 
LUMO orbital for that dimer, see Fig. [T31 By sliding 



TABLE IV. The most important electronic couplings and the 
reorganization energy in the PBI-F2 crystal for electrons and 
holes, cf. Fig. [Til 





h+ [meV] 


e [meV] 


h+ [meV] Ref. 97 e" 


[meV] Ref. 97 




0.251 


129.234 


2 


107 


V2 


2.398 


0.452 






V3 


0.010 


0.017 






Vi 


0.003 


0.004 






Vs 


0.001 


0.002 






A 


213 


303 


215 (213) 


309 (307) 




FIG. 11. (Color online) The most important hopping paths 
in the PBI-F2 crystal. 



one monomer relative to the other along the long axis, 
the coupling for holes oscillates depending on the dis- 
placement around zero (22 because the overlap of the two 
HOMO orbitals with same and different phase alternate. 
All the other coupling constants do not differ significantly 
for the two types of charge carriers. This sole difference 
in the coupling results in a maximum electron mobility 
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FIG. 12. (Color online) The mobility for holes (top) and elec- 
trons (bottom) in the PBI-F2 crystal in all three dimensions. 
The parameters are F = 10^ V/m, T = 300 K. 



that is two orders of magnitude bigger than the maximum 
hole mobility, which is achieved in c direction. However, 
in the plane perpendicular to 6, the hole mobility is two 
orders of magnitude bigger than that of electrons, see 

Fig.m 

The calculated reorganization energies, 303 meV for 
electrons and 213 meV for holes, is bigger than those 
for rubrene and pentacene. The values are in very 
good agreement with reorganization energies calculated 
by Delgado et al^ (309 and 307 meV for electrons, 215 
and 213 meV for holes). 

In order to test our master equation approach, some 
calculations were verified with Monte Carlo calculations. 
The results of both methods agree very well within the 
error bars of the Monte Carlo method. As an example 
Fig. [T5I shows the mobility of PBI-F2 in the ah plane cal- 
culated with both approaches. The Monte Carlo simula- 
tions have run for at least 10 ns and have been averaged 
over at least 100 simulation runs, leading to a relative 
average error of less than 1 %. For this example the mas- 
ter equation approach required about 80.000 times less 
CPU time than the Monte Carlo approach. Thus the 
master equation approach is clearly advantageous as it 
is exact within the numerical accuracy of the computer 
while the Monte Carlo approach contains significant and 



TABLE V. The most important electronic couplings and the 
reorganization energy in the PBI-(C4Fg)2 crystal for elec- 
trons, cf. Fig. [161 





e [meV] 


e" [meV] Ref.-^ 




97.7 


95.7 


V2 


33.7 


35.0 




2.1 


2.2 


Vi 


1.1 


0.9 


X 


339 


360 



slowly converging statistical errors. 

4. PBI-(CaF9)2 

A further fiuorinated perylene bisimide was investi- 
gated which was described by Li et al^ The four most 
important electronic couplings are listed in Tab. |V] and 
depicted in Fig. [161 In contrast to the other molecules 
it is striking that there is no symmetry-caused degener- 
ation of the electronic couplings. It is furthermore im- 
portant to notice that the intra-column couplings Vi and 
V2 along the tt stacks, which are parallel to the a axis, 
differ by a factor of 3. This leads to a "trapping" of 
the charge carrier between the monomers which are cou- 
pled by Vi as described in Sec. Ill Bl After jumping from 
one monomer to the next one along Vi, the charge car- 
rier is more likely to jump back to the first monomer 
than to move on along V2. To illustrate this trapping a 
charge trajectory along the a axis, simulated by Monte 
Carlo, is drawn in Fig. [ITl (top). One clearly sees that 
the charge carrier very often oscillates between two sites 
which lowers the mobility of the charge along the stacks. 
For comparison, a charge trajectory in PBI-F2 along the 
high mobility axis is also depicted. No oscillatory mo- 
tions can be found there. 

This peculiarity of PBI-(C4F9)2 becomes important 
when calculating the mobility: Because of the "trapping" 
that is caused by these oscillations, the mobility calcu- 
lated with Eq. ([T5I) or and the Einstein relation PT|) 
is severely overestimated, see Fig. [151 The green dot- 
ted curve is calculated without external field with the 
master equation along with Eq. (fTSt or respectively, 
which is often used in literature. The red solid curve 
is also obtained by the master equation but the direct 
equation for the mobility, Eq. ([TO)) . was applied. The 
maximum mobility between these two curves differ by a 
factor of 2.4. Besides that, the calculation using the dif- 
fusion coefficient and the Einstein relation even results in 
a wrong angle for the maximum mobility. To prove that 
the result of Eq. PH)) (red solid line) is the right one, 
Monte Carlo simulations were conducted (blue points). 
The simulations ran for 10ns and {x) and {{x — 
were averaged over 1000 trajectories. The relative aver- 
age error was about 0.4 % and the deviation of the master 
equation from Monte Carlo was about 0.2 %. The differ- 
ences in the results of Eq. or p^ and (flUl) are not 
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b 



FIG. 13. (Color online) The mobility for electrons and holes in the PBI-F2 crystal in the ab plane (left), ac plane (middle) and 
be plane (right). The parameters are F = 10^ V/m, T = 300 K. 





FIG. 14. (Color online) The PBI-F2 HOMO (left) and the 
LUMO (right) orbital for the dimer which is built by a b shift 
and leads to the coupling Vi, compare Tab. II VI and Fig. 1111 



caused by the electric field. This is shown by the black 
dashed line which was calculated with Eq. ([T5|) but with 
the same field as for the red solid line. One clearly sees 
that the black dashed line does not coincide with the 
red line but with the green line (calculated without field) 
instead, proving that this approach cannot be applied. 




mobility (i [cm /(Vs)] 

FIG. 15. (Color online) Comparison of master equation and 
Monte Carlo results for the electron mobility in PBI-F2 in the 
ab plane. The parameters are F = lO'^ V/m, T = 300 K. The 
two methods show very good agreement. 
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b viaD{F = 10'V/m) - - 
f direct {F = 10^ V/m) 

Monte Carlo • 




FIG. 16. (Color online) The most important hopping paths 
in the PBI-(C4F9)2 crystal. 




time t [ps] 

FIG. 17. Projection of the charge trajectory onto the re- 
spective direction with the highest mobility for PBI-(C4F9)2 
(a direction, top) and PBI-F2 {b direction, bottom) . The pa- 
rameters are F = lO'' V/m, T = 300 K. 



IV. SUMMARY AND CONCLUSIONS 

A quantum chemical protocol for calculating the 
charge carrier mobilities in organic semiconductor crys- 
tals was presented. A hopping model using Marcus the- 
ory has been implemented by means of the master equa- 
tion approach which is more than four orders of magni- 
tude faster than the Monte Carlo method and free from 
statistical errors. In contrast to the master equation, the 
Monte Carlo approach allows to simulate the transport 
parameters with a time dependent framework. However, 
since this is a stochastic method many simulation runs 
are needed in order to achieve an acceptable statistical 



— / 

















a 



0.2 0.4 0.6 
mobility p. [cm^/(Vs)] 



FIG. 18. (Color online) Comparison of the mobility in the ah 
plane of PBI-(C4F9)2 calculated via the diffusion coefficient 
(Eq. HH) and the Einstein relation (Eq. Hill) for F = 
(green, dotted) and F — 10^ V/m (black, dashed), calculated 
directly (Eq. (I10|) . red, solid) and calculated with Monte Carlo 
(Eq. inil, blue points) for F = lO'^V/m (T = 300 K in all 
cases) . 



error. Furthermore, it is important to make sure that 
the stationary state is obtained within the simulation 
time. This is a serious problem for disordered materials. 
Solving the matrix equation ^ describing the stationary 
state instead by means of analytic numerical methods 
guarantees the stationary solution. 

The mobility is often calculated without external field 
and without the master equation by calculating the diffu- 
sion coefficient and applying the Einstein relation. How- 
ever it can easily happen that the diffusion coefficient is 
overestimated in amorphous materials and even in per- 
fect crystals due to to a "trapping" of the charge between 
energetically similar sites. That is why it is more appro- 
priate to calculate the mobility by means of the master 
equation from the charge drift velocity. The obtained re- 
sults fit perfectly with those of Monte Carlo simulations. 
It is advisable even to calculate the diffusion coefficient 
out of the mobility by applying the Einstein relation, be- 
cause in the Eq. (jlOp for the mobility, the trapping can- 
cels. It was shown that the Einstein relation even holds 
for extremely energetically disordered materials for not 
too high electric fields. 

The angular dependence of the mobility in pentacene, 
rubrene, PBI-F2 and PBI-C4F9 was calculated and the 
results were correlated with the morphology of the crys- 
tals. The results for pentacene and rubrene show a good 
qualitative agreement with experimental data. However, 
the absolute values of the mobilities are strongly over- 
estimated as he assumption of localized charge carriers 
that move in a hopping process without any interaction 
with nonlocal lattice vibrations is not completely ade- 
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quate for organic crystals. Nevertheless this simple model 
allows for qualitative transport property predictions. It 
was shown that PBI-F2 appears to be an almost one di- 
mensional n-type semiconductor. 
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Appendix A 

Equation (jlSp is valid even if an external field is applied 
because in this approach the resulting drift is not caused 
by different jump distances parallel or antiparallel to the 
field respectively since these distances rji are fixed by the 
monomer positions. Instead the field influences the jump 
rates Vji, cf. Eq. ([1]). The drift contribution to the jump 
rate would have to be added to or subtracted from the 
actual rate respectively. However, since Uji influences the 
diffusion linearly, the drift cancels when summing across 
all lattice sites. In order to verify this we have computed 



the solution of the time dependent master equation 

^p = Np, (Al) 

which reads 

p(t) = 5]cV'*, (A2) 

i 

with the eigenvalues U and the respective eigenvectors q. 
The diffusion constant can now be calculated via 

D^~{{x{t)~{x){t)f) 

= E {p^it)x, - Y,pAt)^}j , (A3) 

where the summation is across all sites which positions 
are Xi. We have used the Gaussian disorder model de- 
scribed in Sec. Ill El using the Miller- Abrahams hopping 
rate, Eq. (f^ . for the entries of the matrix N, cf. Eq. ([5|). 
Additionally we have used a simple biased random walk 
where the mobility and the diffusion can even be cal- 
culated analytically. Our calculations confirmed that 
Eq. (fTSi) leads to exactly the same results as Eq. (|A3|) as 
long as there is no energetic disorder, i.e. q{E) — 5{E), 
cf. Eq. (PT|) . The reason for the deviations in the case 
(T 7^ have already been explained in detail in Sec. Ill Bl 
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